home *** CD-ROM | disk | FTP | other *** search
/ PC-Blue - MS DOS Public Domain Library / PC-Blue MS-DOS Public Domain Library - NYACC.iso / vol401 / numdiff.arc / DDBAR.FOR < prev    next >
Encoding:
Text File  |  1986-01-03  |  3.9 KB  |  137 lines
  1. C
  2. C     ..................................................................
  3. C
  4. C     SUBROUTINE DDBAR
  5. C
  6. C     PURPOSE
  7. C     TO COMPUTE, AT A GIVEN POINT X, AN APPROXIMATION Z TO THE
  8. C     DERIVATIVE OF AN ANALYTICALLY GIVEN FUNCTION FCT THAT IS 11-
  9. C     TIMES DIFFERENTIABLE IN A DOMAIN CONTAINING A CLOSED INTERVAL -
  10. C     THE SET OF T BETWEEN X AND X+H (H POSITIVE OR NEGATIVE) - USING
  11. C     FUNCTION VALUES ONLY ON THAT INTERVAL.
  12. C
  13. C      USAGE
  14. C     CALL DDBAR(X,H,IH,FCT,Z,)
  15. C     PARAMETER FCT REQUIRES AN EXTERNAL STATEMENT
  16. C
  17. C     DESCRIPTION OF PARAMETERS
  18. C     X   - THE POINT AT WHICH THE DERIVATIVE IS TO BE COMPUTED
  19. C           X IS IN DOUBLE PRECISION
  20. C     H   - THE NUMBER THAT DEFINES THE CLOSED INTERVAL WHOSE END-
  21. C           POINTS ARE X AND X+H (SEE PURPOSE)
  22. C           H IS IN SINGLE PRECISION
  23. C     IH  - INPUT PARAMETER (SEE REMARKS AND METHOD)
  24. C           IH NON-ZERO - THE SUBROUTINE GENERATES THE INTERNAL
  25. C                 VALUE HH
  26. C           IH    =     0 - THE INTERNAL VALUE HH IS SET TO H
  27. C     FCT - THE NAME OF THE EXTERNAL DOUBLE PRECISION FUNCTION
  28. C           SUBPROGRAM THAT WILL GENERATE THE NECESSARY FUNCTION
  29. C           VALUES.
  30. C     Z   - RESULTING DERIVATIVE VALUE - DOUBLE PRECISION
  31. C
  32. C     REMARKS
  33. C     (1)  IF H = 0, THEN THERE IS NO COMPUTATION.
  34. C     (2)  THE (MAGNITUDE OF THE) INTERNAL VALUE HH, WHICH IS DETER-
  35. C          MINED ACCORDING TO IH, IS THE MAXIMUM STEP-SIZE USED IN
  36. C          THE COMPUTATION OF THE ONE-SIDED DIVIDED DIFFERENCES (SEE
  37. C          METHOD.)    IF IH IS NON-ZERO, THEN THE SUBROUTINE GENERATES
  38. C          HH ACCORDING TO CRITERIA THAT BALANCE ROUND-OFF AND TRUN-
  39. C          CATION ERROR.  HH ALWAYS HAS THE SAME SIGN AS H AND IT IS
  40. C          ALWAYS LESS THAN OR EQUAL TO THE MAGNITUDE OF H IN AB-
  41. C          SOLUTE VALUE, SO THAT ALL COMPUTATION OCCURS IN THE CLOSED
  42. C          INTERVAL DETERMINED BY H.
  43. C
  44. C     SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
  45. C     THE EXTERNAL FUNCTION SUBPROGRAM FCT(T) MUST BE FURNISHED BY
  46. C     THE USER. FCT(T) IS IN DOUBLE PRECISION
  47. C
  48. C     METHOD
  49. C     THE COMPUTATION OF Z IS BASED ON RICHARDSON'S AND ROMBERG'S
  50. C     EXTRAPOLATION METHOD AS APPLIED TO THE SEQUENCE OF ONE-SIDED
  51. C     DIVIDED DIFFERENCES ASSOCIATED WITH THE POINT PAIRS
  52. C     (X,X+(K*HH)/10)K=1,...,10.  (SEE FILLIPI, S. AND ENGELS, H.,
  53. C     ALTES UND NEUES ZUR NUMERISCHEN DIFFERENTIATION, ELECTRONISCHE
  54. C     DATENVERARBEITUNG, ISS. 2 (1966), PP. 57-65.)
  55. C
  56. C     ..................................................................
  57. C
  58.       SUBROUTINE DDBAR(X,H,IH,FCT,Z)
  59. C
  60. C
  61.       DIMENSION AUX(10)
  62.       DOUBLE PRECISION X,FCT,Z,AUX,A,B,C,D,DH,HH
  63. C
  64. C     NO ACTION IN CASE OF ZERO INTERVAL LENGTH
  65.       IF(H)1,17,1
  66. C
  67. C     GENERATE STEPSIZE HH FOR DIVIDED DIFFERENCES
  68.     1 C=ABS(H)
  69.       B=H
  70.       D=X
  71.       D=FCT(D)
  72.       IF(IH)2,9,2
  73.     2 HH=.5D-2
  74.       IF(C-HH)3,4,4
  75.     3 HH=B
  76.     4 HH=DSIGN(HH,B)
  77.       Z=DABS((FCT(X+HH)-D)/HH)
  78.       A=DABS(D)
  79.       HH=1.D-2
  80.       IF(A-1.D0)6,6,5
  81.     5 HH=HH*A
  82.     6 IF(Z-1.D0)8,8,7
  83.     7 HH=HH/Z
  84.     8 IF(HH-C)10,10,9
  85.     9 HH=B
  86.    10 HH=DSIGN(HH,B)
  87. C
  88. C     INITIALIZE DIFFERENTIATION LOOP
  89.       Z=(FCT(X+HH)-D)/HH
  90.       J=10
  91.       JJ=J-1
  92.       AUX(J)=Z
  93.       DH=HH/DFLOAT(J)
  94.       DZ=1.E38
  95. C
  96. C     START DIFFERENTIATION LOOP
  97.    11 J=J-1
  98.       C=J
  99.       HH=C*DH
  100.       AUX(J)=(FCT(X+HH)-D)/HH
  101. C
  102. C     INITIALIZE EXTRAPOLATION LOOP
  103.       D2=1.E38
  104.       B=0.D0
  105.       A=1.D0/C
  106. C
  107. C        START EXTRAPOLATION LOOP
  108.       DO 12 I=J,JJ
  109.       D1=D2
  110.       B=B+A
  111.       HH=(AUX(I)-AUX(I+1))/B
  112.       AUX(I+1)=AUX(I)+HH
  113. C
  114. C        TEST ON OSCILLATING INCREMENTS
  115.       D2=DABS(HH)
  116.       IF(D2-D1)12,13,13
  117.    12 CONTINUE
  118. C        END OF EXTRAPOLATION LOOP
  119. C
  120. C        UPDATE RESULT VALUE Z
  121.       I=JJ+1
  122.       GO TO 14
  123.    13 D2=D1
  124.       JJ=I
  125.    14 IF(D2-DZ)15,16,16
  126.    15 DZ=D2
  127.       Z=AUX(I)
  128.    16 IF(J-1)17,17,11
  129. C        END OF DIFFERENTIATION LOOP
  130. C
  131.    17 RETURN
  132.       END
  133. UE
  134. C        END OF EXTRAPOLATION LOOP
  135. C
  136. C        UPDATE RESULT VALUE Z
  137.       I=JJ+1
  138.